% ====================================================== %
% =================== HF Floating Point Design =================== %
% ====================================================== %
% ====================================================== %


clear
clc

rand('state',0);
randn('state',0);


%%
% ************************************************************* %
% ******************* System Main Parameters *********************** %
% ************************************************************* %
N_Pkts = 50;

Rate = 9600 ;                       % Supported rates are (3200-4800-6400-8000-9600) (Note: Rate 12800 is not supported yet)
Interleaver_L = 72;                % Interleaver sizes are (1-3-9-18-36-72)

ESNO_db = 33; %[28 30 33 35];

NOISE_FLAG = 1;                    % 0 for noise free and 1 for AWGN noise
Fractional_Flag = 1;                % zero for symbol space
SAMPLING_FLAG = 1;                 % set it to one for adding sampling offset effect
with_carr = 1;
InitNoiseExist = 0;
N1_Synch = 25;
N2_Synch = 5;

% Carrier Frequency Offset, max frequency offset in the carrier is 0.018(Hz)
Freq_Offset = 0.0;

% The sampling offset shift between Tx and Rx (This will give 0.24 sampling error)
sampling_offset = -1/100000;
%%
% ************************************************************* %
% ******************* HF_INPUT_PARAMETERS *********************** %
% ************************************************************* %
Oversample = 4;
Symbol_Rate = 2400;
Carrier_Frequency = 1800;
Sampling_Frequency = Symbol_Rate*Oversample;
SamplingFrequencykHz = Sampling_Frequency*1e-3;
% ************************************************************* %

% ************************** CHANNEL *************************** %
EsNo = 10.^(ESNO_db/10);

fadingMethod=1;
% fadingMethod:
% 0 = No fading
% 1 = Uniform Power
% 2 = exponentilal Power profile

pathDelaysMs = [0, 2];

% fadding (two segma) bandwidth
DopplerSpreadHz=1.0;

% factor to compensate the equalizer transfer function at the high SNR
equalizerOffset=0.0;


% for sampling error detection
Delay_Samples = 0;

% Number of symbols of the raised cosine filter.
Filter_Length = 12;

% Rolloff factor that satisfies the standard requirements on the filter (by trail)
Rolloff = 0.385;

% Generate the Raised cosine filter
HF_SRRC = rcosfir(Rolloff,Filter_Length/2,8,1,'sqrt');

% the number of the repeated Initial preamble (for AGC equirements)
% also we may not need to make AGC, but we have to detect
% the repeated sequence at the receive
N_Initial_Preamble = 1;

modulationIndex = find(Rate == [3200 4800 6400 8000 9600])+1;

% Constraint length of the encoder
K = 7;

% Code Rate of the encoder
codeRate = 3/4;
% ************************************************************* %


global LUT_PSK_Q_4 LUT_PSK_Q_8 LUT_QAM_16 LUT_QAM_32 LUT_QAM_64 Sync_Preamble Mini_Probe_Plus Mini_Probe_Minus LUT_S LUT_D LFSR_out Load_Location Mini_Probe_Plus_Modulated trellis_t LUT_Current

load Inputs/Preamble       % by (Preamble_LUTS.m), the output variables are (LUT_D,LUT_S,Mini_Probe_Minus,Mini_Probe_Plus,Sync_Preamble)
load Inputs/LUT_PSK_Q_4    % by (MODULATION_LUTS.m)
load Inputs/LUT_PSK_Q_8    % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_16     % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_32     % by (MODULATION_LUTS.m)
load Inputs/LUT_QAM_64     % by (MODULATION_LUTS.m)
load Inputs/alpha          % by (alpha_generation.m)
load Inputs/trellis_t      % by file (generate_trellis.m)
eval(['load Inputs/' 'LFSR_NSymbols_' num2str(modulationIndex)]);  % by (Gen_LFSRs.m), the variable output will be always (LFSR_out)
eval(['load Inputs/' 'Load_Location_Rate' num2str(Rate) '_Interleaver' num2str(Interleaver_L)]); % by (InterleaverLUTs_Generation.m), the variable output will be always (Load_Location)
LUT_Current = Get_Constellation(Rate);
load Inputs/Initial_Noise
load Inputs/init_noise_indecies
load Inputs/init_noise

disp('************************************************************')
disp('******************** INITIALIZATIONS ***************************')
disp('************************************************************')
disp(' ')
disp(['N_Pkts  =  ' num2str(N_Pkts)])
disp(['Rate  =  ' num2str(Rate)])
disp(['Interleaver_L  =  ' num2str(Interleaver_L)])
disp(['AWGN NOISE  =  ' num2str(NOISE_FLAG)])
disp(['ESNO db  =  ' num2str(ESNO_db)])
disp(['FADING  =  ' num2str(fadingMethod)])
disp(['DopplerSpreadHz  =  ' num2str(DopplerSpreadHz)])
disp(['Freq_Offset  =  ' num2str(Freq_Offset)])
disp(['sampling_offset  =  ' num2str(sampling_offset*SAMPLING_FLAG)])
disp('************************************************************')
disp(' ')
for Pkt = 1:N_Pkts
    
    disp(['Packet Number = ', num2str(Pkt)]);
    disp('============================');
    
    %%
    % ***************** INPUT DATA GENERATION ********************** %
    IN_Frames =72;
    Block_Size_In = HF_Get_Blocksize(Rate,Interleaver_L);
    Input_Bits = randint(1,ceil(Block_Size_In*IN_Frames/Interleaver_L));
    % *********************************************************** %
    
    %%
    % *********************************************************** %
    % ********************** HF_TX_Main *************************** %
    % *********************************************************** %
    % BLOCKING
    % ------------
    [M_OUT_Blocks N_Blocks] = HF_Blocking(Input_Bits,Rate,Interleaver_L);  % put the data into blocks
    Input_Bits_Actual=reshape(M_OUT_Blocks.',1,[]);
    
    % BLOCK PROCESSING
    % ------------------------
    Frame_Count = 1;
    for block_index = 1:N_Blocks
        OUT_Bits_Encoded(block_index,:) = HF_Convenc_Tailbiting_K7_Punct_R43(M_OUT_Blocks(block_index,:),trellis_t);%#ok
        OUT_Bits_Interleaved(Frame_Count:Frame_Count+Interleaver_L-1,:) = HF_Interleaving(OUT_Bits_Encoded(block_index,:),Rate,Interleaver_L);
        Frame_Count = Frame_Count+Interleaver_L;
    end
    N_Frames = Frame_Count-1;
    
    % FRAME PROCESSING
    % -------------------------
    for frame_index = 1:N_Frames
        OUT_Symbols(frame_index,:) = HF_Scrambling( OUT_Bits_Interleaved(frame_index,:),Rate);%#ok
        Modulated_Data(frame_index,:) = HF_Modulate(OUT_Symbols(frame_index,:),Rate);%#ok
    end
    
    % FRAMMINGH
    % -----------------
    TX_OUT = 1*HF_Framming(Modulated_Data,Rate,Interleaver_L,N_Frames,N_Initial_Preamble);
    
    % Oversample
    % --------------
    TX_OUT_Oversampled = upsample(TX_OUT,Oversample);
    
    % SRRC Filtering
    % -----------------
    TX_OUT_Filtered = conv(TX_OUT_Oversampled,HF_SRRC);
    % *********************************************************** %
    TX_OUT_Filtered_Passband = TX_OUT_Filtered;
    %%
    % *********************************************************** %
    % ********************** PassBand *************************** %
    % *********************************************************** %
    % put the baseband signal on a carrier (1800 Hz) as mentioned in the
    % standard, also put the frequency offset (carrier frequency offset)
    
    if with_carr == 1
        Carrier_Phase_Tx = 2*pi*((Carrier_Frequency+Freq_Offset)/Sampling_Frequency);
        TX_OUT_Filtered_L = length(TX_OUT_Filtered);

        TX_OUT_Filtered_Passband = sqrt(2)*real(TX_OUT_Filtered .* exp(1i * Carrier_Phase_Tx * (1 : TX_OUT_Filtered_L)));
    end
    % to normalize the power of the Tx signal (also amplitude) we have to
    % multiply it by sqrt(2) to cancel the power of the sinusoidal signal (carrier)
    % *********************************************************** %
    
    %%
    % *********************************************************** %
    % ********************* CHANNEL EFFECT ************************ %
    % *********************************************************** %
    
    % Fading Effect
    % ---------------
    [ RX_IN, fadingGains ] = WattersonFading( TX_OUT_Filtered_Passband, SamplingFrequencykHz, ...
        DopplerSpreadHz, pathDelaysMs, fadingMethod, Oversample, Fractional_Flag);
    
    
    % Noise
    noise = (randn(1,length(RX_IN))+1i*randn(1,length(RX_IN)))/sqrt(2);
    % *********************************************************** %
    
    %%
    for SNR_Index =1:length(ESNO_db)
        tic
        
        disp(['SNR = ', num2str(ESNO_db(SNR_Index))]);
        disp('----------------------------------------------');
        
        %%
        % ******************* ADD AWGN ***************************** %
        Sigma = NOISE_FLAG / sqrt(EsNo(SNR_Index));
        RX_Faded_Noisy = RX_IN + Sigma/alpha(modulationIndex) * noise ;
        
        if InitNoiseExist == 1
            Current_Indecies = init_noise_indecies(Pkt)+1:init_noise_indecies(Pkt+1);
            RX_Faded_Noisy = [Initial_Noise(Current_Indecies) RX_Faded_Noisy];
        end
        % ********************************************************* %
        if SAMPLING_FLAG == 1
            RX_Faded_Noisy = interp1(1:length(RX_Faded_Noisy), RX_Faded_Noisy,...
                1:(1+sampling_offset):length(RX_Faded_Noisy),'spline');
        end
        % ******************* BASEBAND ***************************** %
        % down conversion for the passband signal to be base band again
        
        if with_carr == 1
            Carrier_Phase_Rx = 2*pi*(Carrier_Frequency/Sampling_Frequency);
            RX_Faded_Noisy_L = length(RX_Faded_Noisy);
            TX_OUT_Filtered_Passband_L = length(TX_OUT_Filtered_Passband);
            
            RX_Faded_Noisy_Baseband =sqrt(2)* (RX_Faded_Noisy .* cos(Carrier_Phase_Rx*(1 : RX_Faded_Noisy_L)) -...
                1i * RX_Faded_Noisy .* sin(Carrier_Phase_Rx * (1:RX_Faded_Noisy_L)) );

            Original_RX_IN_Baseband =sqrt(2)* (TX_OUT_Filtered_Passband .* cos(Carrier_Phase_Rx*(1 : TX_OUT_Filtered_Passband_L)) - ...
                1i * TX_OUT_Filtered_Passband .* sin(Carrier_Phase_Rx * (1:TX_OUT_Filtered_Passband_L)) );
        else
            RX_Faded_Noisy_Baseband = RX_Faded_Noisy;
            Original_RX_IN_Baseband = TX_OUT_Filtered_Passband;
        end
        % ********************************************************* %
        %%

        % Raised Cosine filter at the receiver
        RX_IN_Filtered = conv(RX_Faded_Noisy_Baseband,HF_SRRC);
        
        % Original Tx Symbols passed through the Tx and Rx filters
        Original_RX_IN = conv(Original_RX_IN_Baseband,HF_SRRC);
        
        %%
        % ***************** TIME SYNCHRONIZATION ********************* %
        EsNo_test = (10.^(24/10));  % the for the SNR used in filter coefficients calculations (MMSE equation)
        Original_Start = Filter_Length * Oversample;            % 48
        
        [RX_IN_Filtered, First_Sample, Rate_Rx, Interleaver_Rx] = ...
            HF_TimeSynch(RX_IN_Filtered,Oversample, Original_RX_IN, EsNo_test, Original_Start,...
            N1_Synch, N2_Synch);
        
        % Parameters used to test the correctness of the time synchronization.
        First_Sample_Saved(Pkt, SNR_Index) = First_Sample;
        Rate_Saved(Pkt, SNR_Index) = Rate_Rx;
        Interleaver_Saved(Pkt, SNR_Index) = Interleaver_Rx;
        
        Actual_Start = init_noise(Pkt) + Original_Start;
        
        if(First_Sample == Actual_Start)
            FirstSample_Validate(Pkt, SNR_Index) = 1;
        else
            FirstSample_Validate(Pkt, SNR_Index) = 0;
        end
        
%                                 continue
        % ********************************************************* %
        
        %%
        % ************* OFFSET ESTIAMTION AND EQUALIZATION ************ %
        
        % extracting the signs of all mini-probes (71)
        % -------------------------------------------------
        [Mini_Probes_Sign] = Extract_Miniprobes_Sign_Window(Rate,Interleaver_L);
        Start_Sample = First_Sample + (184 + 31 + 41) * Oversample;
        First_Sample_Index = (First_Sample+1) + (184+103)*Oversample;  % first sample of the first data frame
        First_Sym_Index = 1+ (184) + First_Sample/Oversample +103 ;    % the index of the first symbol of the data
        % ********************************************************* %
        
        %%
        % ******************* Equalizer ****************************** %
        
        % Zero Iteration
        % -----------------
        Hard = 1;                        % one to use hard decision, zero for all training
        N1=40;                            % forward filter coefficients
        N2=5;                              % backward filter coefficients
        EsNo_test = (10.^(24/10));  % the for the SNR used in filter coefficients calculations (MMSE equation)
        
        
        second = 0;
        
        [RX_Input_Frames fadingSymbols MSE_Real MSE OUT_Channel_OUT RX_IN_Filtered gg] = DFE_Equalizer(RX_IN_Filtered,...
            Start_Sample,Oversample,N_Frames,Rate,EsNo_test, Original_RX_IN,Hard,...
            N1,N2, Original_Start,SAMPLING_FLAG,second,Pkt);
        
        [OUT_Bits_Decoded DemodDescrambled_Data N_Blocks] = HF_RECEIVER(RX_Input_Frames,N_Frames,...
            Rate,Interleaver_L,Sigma,fadingSymbols,alpha(modulationIndex),MSE_Real, MSE);
        % ********************************************************* %
        
        Return_Tx_Data = reshape(OUT_Bits_Decoded,1,[]);
        [TX_OUT_Filtered,Input_Bits_Actual,OUT_Bits_Interleaved,N_Frames,TX_OUT] = HF_TRANSMITER(Return_Tx_Data,...
            Rate,Interleaver_L,N_Initial_Preamble,alpha,Oversample,HF_SRRC,modulationIndex);
        
        Original_RX_IN = conv(TX_OUT_Filtered,HF_SRRC);
        
        % First Iteration
        % ----------------
        Hard = 0;            % one to use hard decision, zero for all training
        second = 1;
        
        [RX_Input_Frames fadingSymbols MSE_Real MSE OUT_Channel_OUT RX_IN_Filtered gg] = DFE_Equalizer(RX_IN_Filtered, ...
            Start_Sample,Oversample,N_Frames,Rate,EsNo_test, Original_RX_IN, Hard, N1,N2, Original_Start...
            , SAMPLING_FLAG,second,Pkt);
        
        [OUT_Bits_Decoded DemodDescrambled_Data N_Blocks] = HF_RECEIVER(RX_Input_Frames,N_Frames,...
            Rate,Interleaver_L,Sigma,fadingSymbols,alpha(modulationIndex),MSE_Real, MSE);
        % ********************************************************* %
        
        %%
        % *********************** BER Calculations ********************* %
        
        % CHANNEL BER
        % ----------------
        errorPattern=(sign(DemodDescrambled_Data)+1)/2 ~= OUT_Bits_Interleaved;
        errorPattern=reshape(errorPattern',1,[]);
        channelBitErrors(Pkt,SNR_Index)=sum(errorPattern)/length(errorPattern);
        
        % CODED BER
        % --------------
        RX_OUT_Bits = reshape(OUT_Bits_Decoded.',1,Block_Size_In*N_Blocks);
        biterrorvector = Input_Bits ~= RX_OUT_Bits(1:length(Input_Bits));
        BER(Pkt,SNR_Index) = biterr(Input_Bits,RX_OUT_Bits(1:length(Input_Bits)))/length(Input_Bits);
        
        BER;
        tt =BER% [gg BER]
        mean(BER)
        toc
        % ********************************************************* %
    end
end


% ********************* AVERAGE BER ***************************** %
BER_Avg= mean(BER);
channelBitErrors_Avg = mean(channelBitErrors);

channelBitsEbNo=ESNO_db(1:length(channelBitErrors_Avg))+10*log10(1/modulationIndex);
codedEbNo=ESNO_db(1:length(BER_Avg))+10*log10(1/codeRate/modulationIndex);
SNR = ESNO_db(1:length(BER_Avg))+10*log10(1/Oversample);

BER_Avg
channelBitErrors_Avg

figure
semilogy (ESNO_db(1:length(BER_Avg)),BER_Avg)
grid on

eval(['save ' 'results/BER_Avg' ' BER_Avg']);
eval(['save ' 'results/ESNO_db' ' ESNO_db']);
eval(['save ' 'results/BER' ' BER']);
% ************************************************************* %





